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ABSTRACT 

In this paper we develop methods to use a linear optical model to capture the field dependence of wavefront aberrations 
in a nonlinear optimization-based phase retrieval algorithm for image-based wavefront sensing. The linear optical model 
is generated from a ray trace model of the system and allows the system state to be described in terms of mechanical 
alignment parameters rather than wavefront coefficients. This approach allows joint optimization over images taken at 
different field points and does not require separate convergence of phase retrieval at individual field points. Because the 
algorithm exploits field diversity, multiple defocused images per field point are not required for robustness. Furthermore, 
because it is possible to simultaneously fit images of many stars over the field, it is not necessary to use a fixed defocus 
to achieve adequate signal-to-noise ratio despite having images with high dynamic range. This allows high performance 
wavefront sensing using in- focus science data. We applied this technique in a simulation model based on the Wide Field 
Infrared Survey Telescope (WFIRST) Intermediate Design Reference Mission (IDRM) imager using a linear optical 
model with 25 field points. We demonstrate sub-thousandth-wave wavefront sensing accuracy in the presence of noise 
and moderate undersampling for both monochromatic and polychromatic images using 25 high-SNR target stars. Using 
these high-quality wavefront sensing results, we are able to generate upsampled point-spread functions (PSFs) and use 
them to determine PSF ellipticity to high accuracy in order to reduce the systematic impact of aberrations on the 
accuracy of galactic ellipticity determination for weak-lensing science. 
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1. INTRODUCTION 

Previous work has shown that phase retrieval algorithms for image -based wavefront sensing [1,2, 3, 4] benefit greatly 
from using multiple defocused images as data rather than a single in-focus image. Even if an algorithm must work on a 
single image, a fixed defocus is preferable to purely in-focus images. This is partly because defocused images have 
much less dynamic range than in-focus images. This means that for a saturation-limited exposure, a defocused image can 
accept a longer exposure time (or a brighter source) without saturating and capture more total source photons. Thus, it 
achieves a higher average signal-to-noise ratio than an in-focus image. A second advantage is that the defocused image 
has more visible structure and shape over larger scales in the image. Finally, if focus diversity is employed rather than a 
fixed defocus, the robustness of the wavefront sensing algorithm significantly improves in terms of stagnation 
avoidance, stability and accuracy. A specific example of the benefits of focus diversity is the twin image problem: if the 
system’s pupil is symmetric with respect to reflection through the origin, the phase solution from a single image will 
have a fundamental ambiguity: it may correspond either to a given solution or the same solution reflected through the 
origin and with the phase reversed in sign. 

The advantages of focus diversity and defocused images for image-based wavefront sensing make it highly compelling, 
but requiring focus diversity also limits the systems where it can be applied and the operational scenarios where it can be 
useful. The requirement that most scientific applications of a telescope use in-focus images imposes the first 
fundamental limitation: If you require defocus for wavefront sensing, then wavefront sensing and science use of an 
instrument cannot occur at the same time. Since any time spent on engineering activities is a direct loss from science 
time (unless they can be scheduled while the telescope is otherwise unavailable, such as during a slew), the total time 
that can be allocated to wavefront sensing observations is limited. Even once engineering time is available for wavefront 
sensing, the limited time allocation imposes a trade between the number of defocused images and exposure time per 
image. If the wavefront sensing target stars are dim enough that they are not saturating in the available observing time 
(either due to well-depth or saturation avoidance in sample-up-the-ramp detection), this removes the inherent signal-to- 
noise ratio benefit for defocused images. In addition to the time required, there must be a mechanism to introduce the 



defocus (such as a detector motion or weak lenses in a fdter wheel). If wavefront sensing is required frequently, that 
mechanism must be very robust. Finally even if the engineering time and mechanisms are available, defocusing for 
wavefront sensing is incompatible with science that requires continuous observation of the same target over an extended 
period. 

Because of the costs of defocused and focus-diverse wavefront sensing, we are interested in developing an algorithm that 
can work on collections of in-focus images and still achieve the levels of accuracy and robustness we obtain from focus- 
diverse algorithms. In particular we are interested in wavefront sensing for the imager channel of WFIRST as specified 
in the IDRM report [5,6]. One of the major science applications of the WFIRST imager is a weak gravitational lensing 
survey, which requires making shape (ellipticity) measurements of millions of galaxies over a large area of the sky. 
Because the weak-lensing measurement requires very high accuracy, it is very sensitive to systematic errors introduced 
by aberrations in the imaging system changing the shape of PSF. Because it is impractical to design or build an imaging 
system that is corrected to within the required tolerances over a large field of view, it is instead necessary to know the 
system’s field-dependent point-spread function so that it can be removed from the science data by deconvolution. If the 
system aberrations were static over time, it would be sufficient to measure the system aberrations once during the 
commissioning process and use that as a calibration over the life of the mission. In practice, the system’s PSF will vary 
over time due to thermoelastic effects, as the telescope’s thermal environment changes at different pointing angles, and 
additionally vary over much longer time scales as structures outgas and change shape slightly. Again because the 
scientific performance requirements for the weak-lensing measurement are so tight, these small dynamic errors are a 
concern; it is likely to require an ongoing monitoring process, with either continuous monitoring or periodic monitoring. 
Determining the time scale of variation and the required cadence of wavefront monitoring will require further thermal 
and mechanical monitoring which is not yet available, so for this paper we will assume pessimistically that wavefront 
sensing is required either continuously or with such a frequency that requiring dedicated defocused imagery for 
wavefront sensing would impose an unacceptably onerous overhead on science. Thus we conclude that using the primary 
science detectors for defocused wavefront sensing is not a suitable solution. 

An alternate approach would be to use a separate instrument or channel for defocused wavefront sensing. The IDRM has 
only a single imager and two spectrometers (not suitable for high precision wavefront sensing), so this would require 
adding a separate instrument channel for wavefront sensing, with a commensurate increase in cost and complexity. Even 
if an auxiliary instrument suitable for wavefront sensing were in place, it would face several difficulties. The first is non- 
common path aberrations; if it does not share the same tertiary mirror as the imager or adds additional down-stream 
optics, its aberrations will not be the same as the imager’s and any time dependent motion or deformation in the non- 
common elements of either channel will lead to errors in wavefront sensing, which are likely to exceed the acceptable 
limits for weak lensing. Second, an instrument with a narrow field of view will not be able to capture the full field- 
dependent aberrations that effect the imager; this is acceptable if the time -dependent aberrations are primarily uniform 
over the field of view, with the field-dependent aberrations relatively stable. Because of these issues we will not consider 
wavefront sensing by separate instruments. This leaves either dedicated engineering detectors in the imager (with either 
fixed or variable defocus) that can be used independently of the science detectors or using in-focus science data for 
wavefront sensing. Both are potentially reasonable approaches; in this paper we consider the second because it requires 
the least additional hardware. 

Without access to focus diversity, we are left with two questions. First, how to recover the benefits of (phase) diversity 
in algorithmic robustness? Second, how to achieve high signal-to-noise ratios necessary to reach the required accuracy 
requirements imposed by weak-lensing science? In the case of WFIRST the answer to both of these questions is the 
same: exploit WFIRST’s wide field of view. First, because the imager’s field of view is much larger than its isoplanatic 
patch size, the field dependence of aberrations over the field of view introduces a form of phase diversity: two stars from 
different parts of the field of view will exhibit different aberrations from one another. If we can provide a model for this 
dependence, we can exploit it as a diversity mechanism in our algorithms. Second, again because the field of view is 
large, there are many bright unresolved stars in the field of view. If we can jointly fit all or many of these stars rather 
than a handful of defocused images from a single star, we can achieve very large gains in overall signal-to-noise ratio, 
allowing us to achieve high accuracy in the wavefront fits. 

In order to realize the benefits of field diversity, we need a parameterized model of the field dependence of the system 
aberrations; if each star were allowed to have a hilly independent wavefront realization, there would be no gain in 
robustness or signal-to-noise ratio. Since we are considering a particular system and have access to the optical design, a 
prescription retrieval [7] model seems appropriate. That is, we parameterize the field dependence in terms of the 



mechanical degrees of freedom (both rigid body and figure error) of the optical surfaces in the system, and we search for 
a mechanical configuration that is consistent with the observed images. Because the number of mechanical degrees of 
freedom is far smaller than the number of stars times a reasonable number of low-order wavefront terms, this represents 
a large reduction in search space compared to the independent image problem. 

Having selected the mechanical degrees of freedom as our search space, we must still determine how to model their 
relation to the field dependence. Design and performance assessment for optical systems are typically done with real ray 
trace models of the system. However, non-linear optimization-based phase retrieval algorithms benefit greatly from the 
use of analytic gradient techniques, and it is not obvious how to incorporate a ray trace model into this process without 
imposing a potentially unacceptable computational cost. Because the merit functions in phase diversity are much more 
expensive to compute than a lens design error function, a direct finite difference gradient computation that might achieve 
acceptable performance for the design problem is much more difficult for wavefront sensing. The ideal solution from a 
wavefront-sensing standpoint would be an analytic representation of the field-dependent aberrations of the system, 
produced through an extension of aberration theory to non-rotationally-symmetric and misaligned systems, such as [8] or 
[9, 10]. Such an analysis including both rigid body and figure terms is not available for WFIRST at the time of this 
writing, so we have adopted a linear optical model (LOM) [1 1] to represent the field dependence. 

The LOM is essentially a first-order sensitivity analysis of the optical design; the design is ray traced in an optical design 
code, and then small perturbations are made to the parameters and retraced; this allows the linear dependence of phase in 
the exit pupil of the system to be calculated by finite differences. Note that this is linear in mechanical perturbation, not 
in the functional form of the field dependence of the aberration, which can include higher order dependences on the field 
angle, provided they occur in this small-perturbation limit. Given these sensitivities and the exit pupil aberration maps 
for the configuration space point where the LOM was generated, a field-dependent wavefront error map can be 
generated for a given mechanical realization of the system. For our simulations the LOM was generated at the perfectly 
aligned case, but in practice a priori known static misalignments could be included, giving a LOM that was accurate in 
the region of configuration space near the known point. In order to speed computation in the wavefront sensing process, 
we further projected the first order wavefront sensitivities onto a low order Zernike wavefront basis. For a given field 
point, the linear optical model can be stored as a N x M matrix L n , where N is the number of mechanical degrees of 
freedom (78 in our case) and M is the number of Zernike terms used for the low-order expansion. Thus the forward 
model, converting mechanical configuration points x to a Zernike coefficient map for the nth field point is a simple 
matrix multiplication followed by an addition of the known constant bias (the design residual aberrations, in our case): 

y n - L n x + b n (1) 

If we assume the total joint phase retrieval error metric is given by 

E = Tj w n E n ( 2 ) 

n 

then, given the analytic gradients of E with respect to the vector components of y n , via the chain rule [12] for partial 
derivatives [7A] we can compute the gradient with respect to x by a matrix multiply with the transpose of the linear 
optical model matrix and a sum over the field points; 

y x E = Y jWn L T n (y yn E n ). (2) 

n 

Because the LOM matrix is relatively small compared to the array sizes used to model the optical fields, the addition of 
these two steps in the forward and gradient model for the wavefront sensing algorithm adds negligibly to the total 
computational cost of the algorithm, which is dominated by computing the Fourier transforms needed to propagate 
complex fields between the pupil and image planes. 

The linear optical model formalism provides us a simple and computationally efficient method for multi-field wavefront 
sensing and prescription retrieval. Because it is based on a first-order analysis at a given point in the design space, it is 
not applicable for systems that are expected to see large changes in mechanical state, which would involve higher order 
dependencies on the mechanical configuration not captured in the LOM. However, because the design of the WFIRST 
imager calls for a system that is both well corrected and very stable, we expect that once the telescope has been aligned 



and commissioned on orbit, using the LOM to model small deviations from the static on-orbit alignment should be 
reasonable. 


2. MODEL MATCHING REQUIREMENTS 

Two major factors limit the quality of wavefront retrievals for phase retrieval algorithms: model matching and noise. 
Model matching refers to how well the physical optics simulations in the forward model agree with the actual physical 
world. This depends both on computational and software resources and on knowledge of the system. For example, when 
modeling the physical aperture stop of the system in the computer, it is necessary to discretize the continuous aperture, 
which introduces a mismatch between the model and the real world. This error can be reduced by using shaded edges 
and/or more finely sampling the pupil. Similar issues apply in selecting the sampling in the image domain and in 
discretizing the continuous wavelength band for a broadband system into finite wavelengths. Even given noise-free data, 
these approximations will limit the ultimate accuracy of the algorithm and introduce a trade between desired accuracy 
and computational burden. Because of the combination of high accuracy requirements and the extremely multi-field 
nature of the WFIRST imager, the computational requirements for an adequately model matched phase retrieval model 
for WFIRST will be substantial. For that reason, rather than undertake a comprehensive exploration of the model 
matching requirements, we have considered a single aspect of the model matching problem, discussed below. 

We assessed the impact and sensitivity to model matching by exploring one particular dimension of this problem, the 
sensitivity to the particular choice of wavelengths in a broadband phase retrieval model. Figure 1 shows results of a 
simulation to address this. The vertical axis shows the error in the ellipticity computed from the retrieved point spread 
functions, while the horizontal axis shows the true (total) ellipticity of the point spread functions. Each marker represents 
a single realization of the phase retrieval problem: particular aberrations and noise. We considered monochromatic data 
matched to a monochromatic model, as well as polychromatic models with 10 and 13 discrete wavelengths. For the 
polychromatic wavelengths, we considered both matched and mismatched cases (in the matched cases the same 
wavelengths are used in the data simulation and retrieval; in the mismatched, we simulated with 13 and retrieved with 10 
or simulated with 10 and retrieved with 13). On average the matched cases (Monochromatic data with monochromatic 
retrieval, polychromatic 10 wavelength data with polychromatic 10 wavelength data, and polychromatic 13 wavelength 
data with polychromatic 13 wavelength model) perform similarly well, while all of the mismatched (polychromatic 10 
wavelength data with polychromatic 13 wavelength model and 13 wavelength polychromatic data with 10 wavelength 
polychromatic model ) cases perform comparatively poorly. 

The direct conclusion from this is that 10 discrete wavelengths is not a sufficient to model the point spread functions 
accurately enough that the errors in the phase retrieval results introduced by discretizing the spectra are smaller than the 
noise-limited accuracy at our signal levels: if N wavelengths is sufficient to model the point spread functions with 
required accuracy, than data simulated with more than N wavelengths should be adequately retrieved by a model with 
only N wavelength, which we do not see here. The more general conclusion from this is that is likely that other 
computational aspects of the problem will require high accuracy and that the overall computational difficulty of the 
problem will be significant. 


3. NOISE LIMITED PERFORMANCE 

In the remainder of the paper, we will focus on noise-limited performance. That is, we will ask the question: assuming 
the modeling match issues are all dealt with, at a given noise level, what is the achievable wavefront sensing 
performance? Simulations addressing this topic may use lower fidelity models than would be required to adequately fit 
real-world data by simply making the same approximations in the phase retrieval forward model that are used in the 
simulation of the ‘measured’ data. As an example of this, we adopt throughout a 64x64 grid across the clear aperture of 
the system to reduce computational cost, recognizing that in practice much finer sampling will be required. 

Although for computational reasons we must use simplified models as discussed above, it is important to establish that 
the overall trends we observe in our simulations are consistent with the trends in higher fidelity models. In order to 
investigate this, we have considered simulations using both monochromatic and polychromatic data. We show that in 
both cases, when the phase retrieval model is matched to the data simulation model, the final reconstruction accuracy is 
the same in either case. We also limited the number of field points we considered, both to reduce the computational cost 
of modeling many field points and because only a finite grid of points are computed in our LOM. 



Comparing true PSF ellipticities with ellipti city fits 
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Figure 1. Retrieval accuracy for monochromatic and polychromatic simulations with different numbers of images, [todo: 
revise figures for print] 


However, we would like to able to extrapolate from this scenario to one with a more realistic number of stars. Based on 
photon noise statistics, we expect the residual wavefront error to go as one over the square root of the number of photons 
(for equally bright stars this is proportional to the total number of stars). In order to test this assumption, we simulated 
cases where we used dithering at the existing points to simulate the noise impact of having more than 25 single stars. 
That is, for one field point in our LOM, we generated several (10 in this case, for a total of 250 PSFs) point spread 
functions with different sub-pixel displacements and noise realizations as a proxy for having more than the 25 field 
points computed in our LOM. In all these cases, we observe that the final phase-retrieval accuracy depends on the total 
photon budget in the image plane, not these particular modeling choices, as shown in Figures 2 and 3; this gives us 
confidence that the photon noise limits established in this approximate regime will have relevance to more realistic 
systems. In particular, the performance is consistent with the expected one over square root dependence on the total 
number of photons; in Figure 3 we fit that model to the observed performance points. Figure 4 shows representative 
simulated data for a single field point along with the corresponding phase-retrieval fit, as well as a noise-free version of 
the PSF and upsampled versions of both the true and fit point-spread functions. We can see that despite the aliasing and 
noise in the data for individual field points, in the joint retrieval we have been able to reconstruct the underlying 
diffraction-limited point-spread function with high accuracy. 


4. PHOTON BUDGET CONSIDERATIONS 

The model fit in Figure 3 lets us extrapolate from the cases considered so far to scenarios with larger numbers of stars as 
we would expect in a full flight scenario. This is based on the assumption that all modeling matching requirements are 
met and that photon noise statistics continue to dominate phase retrieval accuracy, so that an accuracy estimate based 
only on the total photon budget is a reasonable estimate. Our model fit for ellipticity accuracy is 1.5 / sqrt(total photons). 
Using a simple star frequency spread sheet for the vicinity of the galactic pole, we can estimate a total of 3x10 s photons 
from stars between 12th and 20th magnitudes in a 120 sec exposure over a 1.6 pm to 2.0 pm band, for an estimated 
ellipticity error of 8.7x10 5 if all available stellar sources are used for wavefront sensing. The IDRM weak-lensing 
program requires knowledge of PSF ellipticity to 4.7x 1 0c 4 , so it is at least plausible that an algorithm of the type 
described in this paper working on a single frame of science data would be able to achieve the required accuracy. This 
should not be taken as a precise prediction of real performance in practice, but rather as a plausibility argument that the 


approach we have outlined in this paper is basically suitable to the problem at hand and that further development of this 
approach and more detailed simulation studies are warranted to establish more realistic error bounds. 


a) 

w 

2 

a: 

& 

g 

ID 

-t: 

a> 


.S' 

o 

Q. 

ED 


Ellipticity versus number of photons 



Monochromatic (25 PSFs) 

r> 

[I] Monochromatic (250 PSFs) 

V 

Polychromatic (25 PSFs) 

* 

- 

$ 


0 1 2 

3 4 5 


T otal number of photons 


8 

XlO 


Figure 2. Retrieval accuracy for monochromatic and polychromatic simulations with different numbers of dithered images 
in terms of total number of photons. The ellipticity uncertainty shown on the vertical axis is the average over a Monte Carlo 
simulation. 
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Figure 3. Retrieval accuracy for monochromatic and polychromatic simulations with different numbers of dithered images, 
fit to inverse square root model we expect based on photon noise. In this figure we have rescaled the horizontal axis so that 
the square root of total number of photons dependence appears a linear. 
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Figure 4. Real and recovered point spread functions. Note that the PSF structure is much more apparent in the upsampled fit 
PSF than the original data, showing the advantage of a wavefront sensing fit over trying to deconvolve with measured PSFs 
directly. All PSFs are shown at intensity to the 0.25 power. 


5. CONCLUSIONS 

We have developed a multi-field phase retrieval algorithm optimized to the requirements of the WF1RST 1DRM imager 
instrument and weak-lensing science program. Our algorithm uses in-focus images of stars over a wide field of view. 
The field-dependent aberrations of the system act as a form of phase diversity and substitute for focus diversity in this 
context. Because of the extremely large field of view for WF1RST, there will be many stars available for wavefront 
sensing in a given exposure and there will be significant change in the wavefront aberrations over the full field of view. 
We show that even with undersampled (for intensity) images and no focus diversity, we can achieve phase retrieval 
accuracy competitive with state of the art focus-diverse phase retrieval results. This is due to large number of 
simultaneously imaged stars available when using the full field of view of a wide-field telescope like WF1RST. We have 
also shown that, based on a simple photon budget calculation, it is plausible that our algorithm can achieve the required 
weak-lensing accuracy working from a single data frame. 
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